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Influence of particle elasticity in shear testers 



o 
o 

(N 
> 

o 

On 



Ph. 



o ■ 

q ; 

O 

^ ■ 

Oh; 



Dirk Kadau^, Dominik Schwesig^'^, Jorg Theuerkauf^ and Dietrich E. Wolf^ * 



Abstract Two dimensional simulations of non-cohesive 
granular matter in a biaxial shear tester are discussed. 
The effect of particle elasticity on the mechanical behavior 
is investigated using two complementary distinct element 
methods (DEM): Soft particle molecular dynamics sim- 
ulations (Particle Flow Code, PFC) for elastic particles 
and contact dynamics simulations (CD) for the limit of 
perfectly rigid particles. As soon as the system dilates to 
form shear bands, it relaxes the elastic strains so that one 
finds the same stresses for rigid respectively elastic parti- 
cles in steady state flow. The principal stresses in steady 
state flow are determined. They are proportional to each 
other, giving rise to an effective macroscopic friction coef- 
ficient which is about 10 % smaller than the microscopic 
friction coefhcient between the grains. 
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Introduction 

In powder technology the material flow properties are usu- 
ally determined experimentally, e.g. using shear testers [Ij. 
The measurements provide input for various phenomeno- 
logical continuum models that have been proposed in or- 
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der to calculate dense granular flows. These models con- 
tain parameters whose connection to the properties of the 
grains is not yet understood. It is the aim of distinct ele- 
ment simulation methods (DEM) to establish the connec- 
tion between the grain scale and the macroscopic behavior 
directly [2115] . 

The stress-strain behavior of a dense granular assembly 
consists of two parts: the rearrangement of the particles 
on the one hand, and their individual elastic or plastic 
deformation on the other. In this paper we address the 
question, how strongly the flow properties are influenced 
by the grain deformations compared to particle rearrange- 
ments. Therefore, two different distinct element methods 
are used: soft particle molecular dynamics modelling elas- 
tic particles (used here: PFC) and contact dynamics (CD) 
to simulate perfectly rigid particles. By comparing the re- 
sults of the two methods the influence of particle elas- 
ticity can be separated from the effect that particle rear- 
rangements have on the macroscopic stress-strain behavior 
found in e.g. the biaxial shear tester considered here. 

We simulate dense granular flow in a biaxial tester, 
which allows for larger displacements than the Jenike 
shear cell. The biaxial shear tester is a rectangular box ^ 
13 in which the material is sheared under constant strain 
rate in one direction and constant stress in a perpendic- 
ular direction while the load plates in the third direction 
are fixed. In this setup one will reach steady state flow 
(constant volume and stress tensor), and the yield locus 
can be determined. 
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Models 

Both models we use simulate the trajectory of each indi- 
vidual particle by integrating Newton's equations. They 
mainly differ in the way, how the contact forces between 
grains are determined. In soft particle molecular dynam- 
ics (PFC) microscopic elastic deformations of each particle 
have to be taken into account: They determine the forces. 
By contrast, in contact dynamics (CD) particles are con- 
sidered as perfectly rigid, and forces are calculated from 
the volume exclusion constraint. 

In order to make the comparison between the two 
methods as stringent as possible, we simulate a very sim- 
ple two dimensional system of noncohesive round particles 
(disks) with Coulomb friction. The boundary conditions 
and initial particle configurations are the same for both 
types of simulations. In both cases the normal restitution 
coefficient is zero and the friction coefficient is /i = 0.3. 
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In order to give specific numbers we also introduce dimen- 
sional parameters, which have no effect on the simulations: 
The average particle radius is 1 mm, the mass density 
of the material the particles are made of is 10^ kg/m'^. 
Whereas no further parameters enter CD, the molecular 
dynamics algorithm requires the specification of two stiff- 
ness parameters as the particles are not perfectly rigid. 



2.1 

Molecular Dynamics Simulation 

The Particle Flow Code used here is based on soft par- 
ticle molecular dynamics 3 ,6 . Deformations of colliding 
particles are represented by the overlap of idealized disks, 
i.e. interaction forces between the particles are functions 
of the overlap (= sum of the two disk radii - distance 
between the disk centers). 

The force between particles during contact is calcu- 
lated with mechanical elements such as springs and dash- 
pots [7||S1 . In Fig. n a basic visco-elastic contact model 
for two particles in contact is depicted schematically. The 
contact force is decomposed into a normal and a tangen- 
tial component. In the simplest case used here the normal 
force is assumed to depend linearly on the overlap (the 
displacement of the spring). The dashpot contributes a 
normal dissipative force proportional to the time deriva- 
tive of the overlap. The sum of both contributions is 
restricted to be repulsive, i.e. tensile normal forces are 
not allowed, because the particles are assumed to be non- 
cohesive. 

The tangential component of the contact force is imple- 
mented in terms of a linear spring-dashpot as well, where 
the displacement of the spring is the integral of the tangen- 
tial relative velocity over the time of nonzero overlap. This 
represents static friction, hence it is restricted to absolute 
values smaller than /j,F„. When this threshold is reached, 
the tangential relative motion is regarded as sliding with 
sliding friction ^F^ (directed opposite to the tangential 
relative velocity). 

The boundary conditions of the system are realized 
by either strain driven walls or stress driven walls without 
any friction. For all of the PFC simulations the commercial 
Particle Flow Code (PFC2d verS.O j2j) was used. 

For the PFC simulations, disks and walls with normal 
and tangential stiffness coefficients of 1 N/m are used. 
This low stiffness has been chosen in order to empha- 
size the possible difference compared to the perfectly rigid 
particles in CD (stiffness oo). Moreover this keeps the 
computing time for PFC low (it increases with the square 
root of the stiffness). The viscous damping coefficient is 
set equal to the critical damping factor at each contact 
so that the restitution coefficients are zero for all particle 
contacts. 
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Contact Dynamics 

The contact dynamics simulation method has been used 
for the simulation of dense and dry granular materials 
since the beginning of the 1980's. Presuming that the be- 
havior of these materials is governed by perfect volume 




Fig. 1. Mechanical contact model used in the PFC simulation: 
Decomposition of contact forces between particles in normal 
(left figure) and tangential (right figure) direction. In both di- 
rections a spring and a dashpot is applied. Tensile forces are 
prevented in normal direction (left figure), in tangential direc- 
tion a frictional slider is used for realizing Coulomb friction 
(right figure). 



normal force 



tangential force * 



Fig. 2. Contact laws in Contact Dynamics: Perfect volume 
exclusion (Signorini graph, left figure) determines the normal 
force _Fjv, the tangential force Ft is determined by exact imple- 
mentation of Coulomb's law of friction with threshold -Ft, max 
(right figure). 



exclusion and static Coulomb friction 
CD implements them as constraints without modelling the 
microscopic elastic grain deformations underlying them. 
Although CD can deal with more general systems as well 
[TB| we restrict ourselves to the simplest case here. 

Volume exclusion of perfectly rigid particles is char- 
acterized by the Signorini graph. Fig. [2Ileft): When two 
particles touch each other {d = 0) the normal force F^ 
at the contact assumes whatever positive value is needed 
to avoid interpenetration. This value can in principle be 
arbitrarily large (Fig. |21 left). Otherwise the particles do 
not interact. The tangential force Ft at a contact is de- 
termined by the Coulomb friction law. Fig. [3 (right): A 
sticking contact (with relative tangential velocity vt = 0) 
can bear any tangential force with absolute value up to 
a threshold -Ft^max = f^Fn- Exceeding this threshold leads 
to sliding of the two particles in contact. If sliding occurs 
(wt ^ 0) the magnitude of the tangential force is Ft, max, 
while its direction is opposite to the sliding velocity. 

For dense systems consisting of many particles with 
a complex contact network, the calculation of all contact 
forces is a global problem, because every force at a contact 
influences all other contact forces. There exists no analyt- 
ical solution for this complex problem |17II18| . so that an 
iterative procedure is applied: The forces at every contact 
are calculated in a random order repeatedly until they all 
comply with each other. In the simulations presented here 
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a fixed number of iterations is chosen, which according to 
our previous experience is large enough that the forces 
have converged within a small tolerance. If the number of 
iterations would be too small, quasielastic effects would 
occur [T^. 

3 

Simulated System 

We simulate a rectangular system confined by friction- 
less walls perpendicular to the x- and y-axes. Fig. U] The 
boundary conditions are such that the material is sheared 
under constant stress a^x at the yielding walls (perpen- 
dicular to the X-axis) and with a constant strain rate iyy 
at the pushing walls (perpendicular to the y-axis). The 
stress at the pushing walls, ayy, and the strain rate at the 
yielding walls, ixx, are evaluated. Since there is no wall 
friction in our simulations, ayy and ctxx are the principal 
stresses. 

In the PFC and CD simulations the same initial config- 
urations are used. The system consists of about 800 round 
particles (disks) with a Gaussian distribution of diameters 
cut off below 1/2 and above 3/2 of the average value. The 
width of the distribution is 0.1 times the average particle 
diameter. The polydispersity is important to avoid layer- 
ing within the system. 

The initial configuration is prepared in the following 
way: The particles are randomly distributed at low con- 
centration at first. Using contact dynamics the system is 
then biaxially compacted under small constant forces on 
the four walls (leading to stresses small compared to the 
ones apphed later). This leads to a dense packing (Fig. 
13 left). The forces in cc-direction were twice as large as 
those in y-direction in order to obtain an aspect ratio of 
about 2:1. If the average particle radius is taken as 1 mm, 
the system size is about 43mm (a;-direction) times 73mm 
(j/-direction) . Whereas in CD this configurations can be 
directly used as initial configuration {t — 0), a further 
preprocessing step is needed for PFC. Using molecular 
dynamics the configuration is compressed a little in x- 
direction until the desired stress axx is reached due to elas- 
tic response of the grains (here: small overlaps), while the 
extent in y-direction is kept constant. Therefore the ini- 
tial volume of the PFC-simulation (Fig.O right) is about 
6% smaller in the x-direction than the one of the CD- 
simulation, but otherwise the initial configurations are the 
same. 

To reduce fluctuations the results were averaged over 
10 similar systems where particles are randomly dispersed 
and located while keeping the same average distribution. 
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Results 

These dense systems (Fig. are sheared under a con- 
stant strain rate iyy = —0.02 s~^ and different values of 
constant stress axx- The stresses are calculated by divid- 
ing the force on the corresponding wall by the wall area, 
where the average disk radius, 1 mm, is taken as size per- 
pendicular to the xy-plane. 

For axx = 10 Pa an example of a system after shear- 
ing for i = 42 s is shown in Fig. ^ As will be discussed 




Fig. 3. Initial configurations used for CD (left) and PFC 
(right) for shearing under fixed strain rate in j/-direction and 
fixed stress in a;-direction. In the right figure the left and right 
walls have been moved together leading to elastic response of 
the particles in order to fix the stress a^x- 




Fig. 4. Typical configuration in steady state flow at the same 
time for CD (upper) and for PFC (lower panel). Both simula- 
tions started with the same configuration (up to a horizontal 
compression by about 6% for PFC, see text). 



below, at this time steady state flow has already been 
established. The vertical size of the simulation boxes is 
the same, because the strain rate iyy is the same in both 
simulations. Remarkably, the horizontal size difference is 
still about 6%, as for the initial configurations, although 
it could have evolved differently depending on the elastic 
properties, and in fact transiently does so. The strain built 
up temporarily due to elasticity vanishes again in steady 
state flow. The final particle positions are different in the 
two simulation methods, although they were initially the 
same, which is no surprise in view of the highly nonlinear 
dynamics. In particular the force network evolves differ- 
ently. On the other hand, Fig.Elcompares the distributions 
of normal forces for a snapshot like Fig. 0| The distribu- 
tions obtained with CD and PFC are similar, thus the 
force statistics in the systems are comparable. 

We evaluated the velocities of the grains for the CD- 
simulations and found shear localization in shear bands. 
In contrast to previous simulations |20II21II22| . where the 
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Fig. 5. Comparison of tlie distribution of normal forces from 
a snapshot during steady state flow for CD and PFC. 



left wall was fixed, i.e. strain controlled, whereas the right 
one was stress controlled, the shear band which forms ini- 
tially, is oriented along either diagonal, breaking the sym- 
metry of our system spontaneously. Later on most of the 
time four shearbands exist which form roughly a parallelo- 
gram around the center which hardly moves. The dynam- 
ics of the shear bands is very complex showing eddies and 
strong fluctuations. This will be analyzed in more detail 
in a forthcoming publication |23j . 

The measured strain rate exx^ averaged over 10 simi- 
lar systems is shown in Fig. |B| One can distinguish three 
time intervals: During the first 2.5 s the elastic particles 
are compressed, i.e. kxx < —^yy = 0.02s~-^. As expected, 
this region is much smaller for CD, where the system of 
rigid disks gets compacted for only about 1 s, mostly due 
to numerical errors creating tiny overlaps, which for per- 
fect convergence of the force iterations would not occur. 
However, a small amount of compaction due to particle 
rearrangements is also possible. The second time inter- 
val extends up to 10 s in both simulation models and is 
characterized by ixx > ~^yy: which means that the sys- 
tem dilates. The volume increases to allow for shearing 
j24| . As already seen in Fig. 01 the elastic system dilates 
more, compensating for the initial compression, so that 
finally the volumes of both systems differ by roughly the 
same percentage as initially. Obviously the elastic energy 
stored during the compression phase is completely released 
during dilation. After this transient the system reaches a 
steady state where the strain rate fluctuates around the 
average value (ixx) = —^yy In this region the average 
volume remains constant. The PFC-results agree quali- 
tatively with the ones obtained with molecular dynamics 
simulations by Luding pnil21ll22j . 

The time dependence of the stress <Tyy is shown in 
Fig- 13 For elastic particles (PFC) it reaches a pronounced 
maximum at the end of the compression phase after about 
3 s. During the dilatancy phase it decreases and finally 
fiuctuates around a constant value in steady state fiow. 
For rigid particles (CD) the stress maximum is not very 
pronounced: In order to distinguish it convincingly from 
the fiuctuations, one would have to average over more than 
10 independent runs. This confirms that the stress maxi- 
mum is due to elastic compression of the system. On the 
other hand, the average steady state values of ayy are 
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Fig. 6. Time dependence of the the strain rate e^x at fixed 
eyy = —0.02 s~^ and a^x = 10 Pa in comparison for PFC (left) 
and CD (right). Strain curves for different fixed axx show the 
same behavior (not shown). 




Fig. 7. Time dependence of the the stress ayy at fixed eyy — 
-0.02 and axx = 10 Pa for PFC and CD. 



the same for PFC and CD within the error bars. Elas- 
tic contributions to ayy cannot be observed any more in 
agreement with our picture, that the elastic deformation 
relaxes during dilation at the onset of shearing. 

We repeated the simulations for four different stresses 
ffxx between 5 and 25 Pa keeping the strain rate eyy = 
— 0.02s~^ the same. The time-dependence of the strain 
rate ixx as well as of the stress ratio Cyy/axx is in all 
cases the same as in Fig. respectively Fig.0 In particu- 
lar, {fJyy) /uxx time-averaged in steady state fiow does not 
depend on axx, see Fig. |H1 ^ 

As both principal values of the stress tensor are pro- 
portional to each other in steady state flow, this implies 
a linear effective yield locus. The effective friction angle </> 
is given by 

, i^yy) ^xx 



{ayy) + ax 



(1) 



With (ayy) /axx ~ 1-71 ± 0.05 this implies an effective 
friction coefficient 

Aicff = tan 0w 0.27 ±0.02 (2) 

which is slightly smaller than the friction coefficient be- 
tween the grains, fi = 0.3. 

The effect of system size was investigated by increasing 
the number of particles to 10000. However, no effect of the 
system size on the stress ratio was found. 

^ Although the average values for CD seem to be slightly 
larger than those for PFC by 1 to 3 %, the error bars are so 
big that the above statement that elasticity does not play a 
significant role in steady state flow is still valid. 
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Fig. 8. The average stress ayy in steady state flow is propor- 
tional to the appUed stress axx- 



Conclusion 

The results presented above for a two dimensional model 
of a biaxial shear test show that particle elasticity only af- 
fects the stress and volumetric strain during a short tran- 
sient. Dilatancy leads to elastic strain relaxation so that 
the stress in steady state flow is essentially determined by 
rearrangements of the particles. This could be proven by 
comparing the simulation results obtained for elastic and 
perfectly rigid particles, respectively. 

In steady state flow the principal stresses turned out 
to be proportional to each other. Using Mohr-Coulomb 
theory we determined the effective macroscopic friction 
coefficient of the granular material, which is the ratio of 
the shear stress to the normal stress at a shear plane. 
We could relate it to the microscopic friction coefflcient 
between the grains, which is slightly larger. 

We also checked the influence of the strain rate iyy 
and found essentially no effect. This can be explained 
by comparing the orders of magnitude of inertia forces 
due to the prescribed strain rate to the forces in the sys- 
tem caused by the applied stress Uxx- The typical time 
scale for inertia effects is given by the inverse strain rate 
i/iyy, the length scale is given by the system size Ly, 
the typical mass is the average particle mass prj] (we ne- 
glect the factor tt here). Thus, typical inertia forces are: 
-Pinertia ~ P't'o-^yi^y Typical forccs duc to the applied 
stress axx are estimated by multiplying with the system 



size Ly (and its thickness ro): Fc: 



(JxxLyra- The in- 



teresting ratio: ^inertia/-P'ext ~ PTQ^^ly/crxx- Inserting the 
used values results in: Jlnertia/-Pext ~ 4 • 10~^ which is a 
small number so that we are well in the region of slow 
(quasistatic) deformations where we do not expect an in- 
fluence of the applied strain rate on the results as is fomid 
here. This argument shows, that only the dimensionless 
ratio pr^^yjoxx matters, as long as elastic effects are neg- 
ligible: Increasing is equivalent to decreasing Oxx- 



References 

1. J. Schwedes. Review on testers for measuring of properties 
of bulk solids. Granular Matter, 5:1-43, 2003. 



2. P.A. Cundall. A computer model for simulating progres- 
sive large scale movements of block rock systems. In Pro- 
ceedings of the Symposium of the International Society of 
Rock Mechanics, volume 1, pages 132-150, Nancy, Prance, 
1971. 

3. P.A. Cundall and O.D.L Strack. A discrete model for gran- 
ular assemblies. Geotechmque, 29:47-65, 1979. 

4. M. Nowak. Spannungs-/Dehnungsverhalten von Kalkstein 
in der Zweiaxialbox. PhD thesis, TU Braunschweig, 1994. 

5. R. J. M. Janssen and H. Zetzener. Measurements on co- 
hesive powder with two biaxial shear testers. Chemical 
Engineering & Technology, 26(2):147-151, ?? 

6. P.A. Cundall. Distinct element models of rock and soil 
structure. In Analytical and computational methods in 
engineering and rock mechanics, London, England, 1987. 
Allen & Unwin. 

7. D.E. Wolf Modelling and computer simulation of gran- 
ular media. In K.H. Hoffmann and M. Schreiber, edi- 
tors. Computational Physics: Selected Methods - Simple 
Exercises - Serious Applications, pages 64-94, Heidelberg, 
1996. Springer. 

8. S. Luding. Molecular dynamics simulations of granular 

materials. In H. Hinricliseri and D. E. Wolf, editors. The 
Physics of Granular Media, Berlin, Germany, 2004. Wiley- 
VCH. 

9. Itasca consulting, www.itasca.com. Particle Flow Code in 
2 Dimensions. Online manual PFC2d version 3.0. 

10. P. Lotstedt. Mechanical systems of rigid bodies subject to 
unilateral constraints. SIAM J. Appl. Math., 42:281-296, 
1982. 

11. M. Jean and J. J. Moreau. Unilaterality and dry friction 
in the dynamics of rigid body collections. In Proceedings 
of Contact Mechanics International Symposium, pages 31- 
48, Lausanne, Switzerland, 1992. Presses Polytechniques et 
Universitaires Romandes. 

12. M. Jean. The non-smooth contact dynamics method. 
Comput. Methods Appl. Engrg., 177:235 257, 1999. 

13. J.J. Moreau. Some numerical methods in multibody dy- 
namics: application to granular materials. Eur J Mech, 
A/Solids, 13(4):93-114, 1994. 

14. T. Unger and J. Kertesz. The contact dynamics method 
for granular media. In Modeling of Complex Systems, pages 
116-138, Melville, New York, 2003. American Institute of 
Physics, cond-mat/0211696. 

15. L. Brendel, T. Unger, and D. E. Wolf Contact dynamics 
for beginners. In H. Hinrichsen and D. E. Wolf, editors. 
The Physics of Granular Media, Berlin, Germany, 2004. 
Wiley- VCH. 

16. D. Kadau, G. Bartels, L. Brendel, and D. E. Wolf Pore 
stabilization in cohesive granular systems. Phase Trans., 
76(4-5):315-331, 2003. 

17. F. Radjai, L. Brendel, and S. Roux. Nonsmoothness, inde- 
terminacy, and friction in two dimensional arrays of rigid 
particles. Phys. Rev. E, 54(1):861, 1996. 

18. T. Unger, J. Kertesz, and D. E. Wolf Force indetermi- 
nacy in the jammed state of hard disks. Phys. Rev. Lett, 
94:178001, 2005. 

19. T. Unger, L. Brendel, D. E. Wolf, and J. Kertesz. Elastic 
behavior in contact dynamics of rigid particles. Phys. Rev. 
E, 65(6):061305, 2002. 

20. S. Luding, R. Tykhoniuk, J. Tomas, L. Heim, M. Kappl, 
and H.-J. Butt. Flow behavior of cohesive and frictional 
fine powders. In Y. Shimizu, R. D. Hart, and P. A. Cun- 
dall, editors. Numerical Modeling in Micromechanics via 
Particle Methods - 2004, pages 157-163. A. A. Balkema, 
2004. PFC Symposium procedings. 



6 



21. S. Luding. Anisotropy in cohesive, frictional granular me- 
dia. J. Phys.: Condens. Matter, 17:S2623-S2640, 2005. 

22. S. Luding. Shear flow modeling of cohesive and frictional 
fine powder. Powder Technology, 2005. In press, corrected 
proof available online 23 May 2005. 

23. J. Torok, D. Kadau, and D. E. Wolf. Properties of shear 
bands in biaxial tests, in preparation. 

24. O. Reynolds. On the dilatancy of media composed of rigid 
particles in contact. Philos. Mag., Ser. 5(20):469, 1885. 



